\section{用SageMath计算举例}

\begin{frame}[fragile]
  \frametitle{基本算术}

  这里我们介绍如何用SageMath来做本章中的运算。SageMath的官网是 \url{https://www.sagemath.org/}.
  MacOS和Linux系统安装很方便，
Windows系统需要通过wsl安装Ubuntu系统后安装。官网不仅有各组件的文档，而且提供了一些入门的参考书。
我们先介绍点基本的算术。

  ~

四则运算：
\begin{minted}[texcl]{sage}
sage: 2+3  # 求和 $2+3$
5
sage: 2-3  # 求差 $2-3$
-1
sage: 2*3  # 乘积 $2\times 3$
6
sage: 2/3  # 求商 $2/3$ 
2/3
sage: 2/3+5/3+3/2  # 求和 $2/3+5/3+3/2$
23/6
sage: 7//3  # 求部分商， $7$除以$3$的整数部分
2
sage: 7%3  # 求 $7$除以$3$的余数
1
sage: divmod(7,3)  # 带余除法（同时求部分商和余数）
(2, 1)
\end{minted}


\end{frame}

\begin{frame}[fragile]
指数对数运算：
\begin{minted}[texcl]{sage}
sage: 2^3  # 求 $2^3$
8
sage: 2**3  # 求 $2^3$
8
sage: 2^(1/2)  # 定义 $2^{1/2}$
sqrt(2)
sage: sqrt(-2)  # 定义 $\sqrt{-2}$
sqrt(-2)
sage: sqrt(-2)^3+1  # 求 $\sqrt{-2}^3+1$ 
-2*sqrt(-2) + 1
sage: (2.0)^(1/2)  # 求$2^{1/2}$的近似值
1.41421356237310
sage: numerical_approx(sqrt(2))  # 求$2^{1/2}$的近似值
1.41421356237310
sage: 2^(1/3)  # 定义 $2^{1/3}$
2^(1/3)
sage: log(5,2)  # 定义$\log_2 5$
log(5)/log(2)
sage: log(5.0,2)  # 求 $\log_2 5$ 的近似值
2.32192809488736
sage: numerical_approx(log(5,2))  # 求 $\log_2 5$ 的近似值
2.32192809488736
\end{minted}

\end{frame}

\begin{frame}[fragile]
整数的$m$-进制展开：

\begin{minted}[texcl]{sage}
sage: 10.digits(base=2)  # $10$的二进制展开：$10 = 2^1 + 2^3$
[0, 1, 0, 1]
sage: 121.digits(base=5)  # $121$的二进制展开：$121 = 1 + 4\times 5 + 4\times 5^2$
[1, 4, 4]
\end{minted}

~

更多：

\begin{minted}[texcl]{sage}
sage: abs(-2)  # 求绝对值
2
sage: floor(sqrt(101))  # 求不大于$\sqrt{101}$ 的整数
10
sage: ceil(sqrt(101))  # 求不小于$\sqrt{101}$的整数
11
sage: factorial(10)  # 求阶乘 $10!$
3628800
sage: 10.binomial(2)  # 求组合数 $\binom{10}{2}$
45
\end{minted}

\end{frame}

\begin{frame}[fragile]
  \frametitle{因子}

Sage提供了方法分解整数，求最大公因子、最小公倍数，求整数所有的因子 (\mintinline{sage}{divisors(n)})，因子的个数 (\mintinline{sage}{number_of_divisors(n)}),
$n$的因子的$k$-次幂之和 (\mintinline{sage}{sigma(n,k)}):

\begin{minted}[texcl]{sage}
sage: 100.factor()  # 因子分解
2^2 * 5^2
sage: divisors(100)  # 求所有因子
[1, 2, 4, 5, 10, 20, 25, 50, 100]
sage: number_of_divisors(100)  # 求因子个数
9
sage: sigma(100,0);sigma(100,2)
9
13671
sage: 100.gcd(55)  # 求 $(100,55)$
5
sage: 100.gcd(55).gcd(15)  # 求 $(100,55,15)$
5
sage: 31.xgcd(13)  # 求 $(31,13)$ 及 Bézout 等式
(1, -5, 12)
sage: 100.lcm(55);100.lcm(55).lcm(30)  # 求 $[100,55], [100,55,30]$
1100
3300
\end{minted}

\end{frame}

\begin{frame}[fragile]
  \frametitle{素数}

Sage中有 \mintinline{sage}{Primes()} 或 \mintinline{sage}{is_prime(n)} 用于素性检验：

\begin{minted}[texcl]{sage}
sage: 2^(2^4)+1 in Primes()
True
sage: 2^(2^5)+1 in Primes()
False
sage: factor(2^(2^5)+1)
641 * 6700417
sage: is_prime(2^(64)+1)
False
sage: all(is_prime(2^(2^n)+1) for n in range(0,5))
True
sage: any(is_prime(2^(2^n)+1) for n in range(5,10))
False
\end{minted}

~

列举素数：

\begin{minted}[texcl]{sage}
sage: primes_first_n(10)  # 前 $10$ 个素数
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
sage: list(primes(50,100))  # $50$到$100$之间的素数
[53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
\end{minted}
\end{frame}

\begin{frame}[fragile]
求某个数之后或之前的素数：

\begin{minted}[texcl]{sage}
sage: next_prime(100)  # $100$之后的第一个素数
101
sage: previous_prime(100)  # $100$之前的第一个素数
97
\end{minted}

~

计算$\pi(x)$:

\begin{minted}[texcl]{sage}
sage: prime_pi(1000.1)  # 不超过 $1000.1$ 的正素数的个数
168
\end{minted}

~

求$p$-进（指数）赋值：

\begin{minted}[texcl]{sage}
sage: v = QQ.valuation(2)
sage: v(1000)  # 求 $1000$ 的 $2$-进赋值
3
sage: v = QQ.valuation(5)
sage: v(1000)  # 求 $1000$ 的 $5$-进赋值
3
\end{minted}
\end{frame}

\begin{frame}[fragile]
  \frametitle{求解Diophantine方程}

可以调用 \mintinline{sage}|solve_diophantine()| 方法：
 
\begin{minted}[texcl]{sage}
sage: R.<x,y,z> = PolynomialRing(ZZ)
sage: solve_diophantine(434*x+182*y-42)  
(13*t_0 - 15, -31*t_0 + 36)
sage: solve_diophantine(434*x+182*y+10*z-4)
(t_0, -212*t_0 + 5*t_1 + 2, 3815*t_0 - 91*t_1 - 36)
sage: sol = solve_diophantine(x^2-y == 0); sol
(t, t^2)
sage: [(sol[0].subs(t=t),sol[1].subs(t=t)) for t in range(-3,4)]
[(-3, 9), (-2, 4), (-1, 1), (0, 0), (1, 1), (2, 4), (3, 9)]
sage: sol = solve_diophantine(x^2 + y^2 == z^2); sol
(2*p*q, p^2 - q^2, p^2 + q^2)
sage: sorted(solve_diophantine(x^2-2*y^2-1), key=str)
[(-sqrt(2)*(2*sqrt(2) + 3)^t + sqrt(2)*(-2*sqrt(2) + 3)^t - 3/2*(2*sqrt(2) + 3)^t - 3/2*(-2*sqrt(2) + 3)^t,
  3/4*sqrt(2)*(2*sqrt(2) + 3)^t - 3/4*sqrt(2)*(-2*sqrt(2) + 3)^t + (2*sqrt(2) + 3)^t + (-2*sqrt(2) + 3)^t),
 (sqrt(2)*(2*sqrt(2) + 3)^t - sqrt(2)*(-2*sqrt(2) + 3)^t + 3/2*(2*sqrt(2) + 3)^t + 3/2*(-2*sqrt(2) + 3)^t,
  -3/4*sqrt(2)*(2*sqrt(2) + 3)^t + 3/4*sqrt(2)*(-2*sqrt(2) + 3)^t - (2*sqrt(2) + 3)^t - (-2*sqrt(2) + 3)^t)]
\end{minted}

可输入 \mintinline{sage}|?Expression.solve_diophantine| 来了解。
\end{frame}

\begin{frame}[fragile]
  \frametitle{连分数}

  \mintinline{sage}|continued_fraction()| 和
  \mintinline{sage}|continued_fraction_list()|
  可用于计算实数的连分数：

  \begin{minted}[texcl]{sage}
sage: (13/27).continued_fraction()  # $13/27$ 的连分数展开
[0; 2, 13]
sage: 0 + 1/(2 + 1/13)
13/27
sage: continued_fraction(22/45)
[0; 2, 22]
sage: continued_fraction(pi)  # 圆周率的连分数展开
[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]
sage: continued_fraction_list(pi, nterms=20)  # 圆周率的连分数展开至20项
[3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2]
sage: continued_fraction(5^(1/3))
[1; 1, 2, 2, 4, 3, 3, 1, 5, 1, 1, 4, 10, 17, 1, 14, 1, 1, 3052, 1, ...]
  \end{minted}
\end{frame}

\begin{frame}[fragile]


  求连分数的第$n$个渐进分数$\frac{p_n}{q_n}$以及$p_n, q_n$:
  \begin{minted}[texcl]{sage}
sage: cf=continued_fraction(pi)
sage: cf.convergent(0)  # 第 $0$ 个渐近分数（即整数部分）
3
sage: cf.convergent(1)  # 第 $1$ 个渐近分数
22/7
sage: cf.convergent(2);cf.convergent(3);cf.convergent(4)
333/106
355/113
103993/33102
sage: cf.p(2);cf.q(2)
333
106
  \end{minted}

  形如$a+b\sqrt{d}$的连分数展开的周期：

  \begin{minted}[texcl]{sage}
sage: K.<sqrt3> = QuadraticField(3)  # 定义二次数域 $\bQ(\sqrt{3})$
sage: continued_fraction(sqrt3).period()  # 求 $\sqrt{3}$ 的连分数展开的周期部分
(1, 2)
sage: cf=continued_fraction(1+2*sqrt3)
sage: cf
[4; (2, 6)*]
sage: cf.period()
(2, 6)
sage: cf.period_length()  # 求 $1+2\sqrt{3}$ 的连分数展开的周期的长度
2
  \end{minted}
\end{frame}

\begin{frame}[fragile]
  \frametitle{编程举例}

下面的函数实现了按照 \S1.3 定理4 的公式
\[
v_{p}(n!)=\left\lfloor\frac{n}{p}\right\rfloor+\left\lfloor\frac{n}{p^{2}}\right\rfloor+\left\lfloor\frac{n}{p^{3}}\right\rfloor+\cdots
\]
求 $v_p(n!)$:

\begin{minted}[texcl]{sage}
sage: def vpfactorial(n,p):
....:     return sum(floor(n/p^k) for k in range(1,floor(log(n,p))+1))
....: 
sage: print(vpfactorial(11,2))
8
\end{minted}

可用直接的计算检验：
\begin{minted}[texcl]{sage}
sage: v = QQ.valuation(2)
sage: v(factorial(11))
8
\end{minted}


\end{frame}
